---
title: "Reproducible_Code_Incentive_Experiment"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(warning = FALSE, message = FALSE)
```

# Setup

```{r, include=TRUE}
library("tidyverse")
library("stargazer")
library("knitr")
library("readxl")
library("tinytex")
library("glmmTMB")
library("margins")
library("gridExtra")
library("xtable")
library("kableExtra")
library("boot")
library("ggsignif")
library("performance")
library("emmeans")
library("lme4")
```

Full dataset (basis for analyses on consent)
```{r}
gsurf_consentexp_analysis_full <- read.csv("Data_Incentive_Experiment.csv") %>% 
  mutate(con_inc_backward = as.factor(con_inc_backward))
```

# Data wrangling
Prepare backward difference coding:
```{r}
#backward difference for factor variable with 4 levels
my.backward.diff <- matrix(c(-3/4, 1/4, 1/4, 1/4, -1/2, -1/2, 1/2, 1/2, -1/4, -1/4, -1/4, 3/4), ncol = 3)

#code contrasts using the matrix defined before for the conditional incentive
contrasts(gsurf_consentexp_analysis_full$con_inc_backward) <- my.backward.diff 
```

Dataset of those who gave consent (basis for analyses on installation)
```{r}
gsurf_consentexp_analysis_full_consent <- gsurf_consentexp_analysis_full %>% 
  filter(consent == 1)

nrow(gsurf_consentexp_analysis_full_consent)
```

Dataset of those who installed (basis for analyses on tracked days)
```{r}
gsurf_consentexp_analysis_full_active <- gsurf_consentexp_analysis_full %>% 
  filter(active_days_threshold >= 1)
nrow(gsurf_consentexp_analysis_full_active)
```

Footnote 3:
```{r}
gsurf_consentexp_analysis_full %>% 
  group_by(study_tracking_num, rec_arm) %>% 
  summarise(installation_rate = mean(installation, na.rm = T)) %>% 
  ungroup() %>% 
  mutate(recruitment = c("ALLBUS - wave 2", "Meta - wave 2", "Meta - wave 1")) %>% 
  select(recruitment, installation_rate)
```

# Descriptive statistics

## For table A2
```{r}
variables <- colnames(gsurf_consentexp_analysis_full %>%
                        select(consent,installation,active_days_threshold, sex, age, educ, german,married,hhsiz,
                               pa_index, kh_index))

table_descriptives1 <- data.frame(Variable = c("Consent", "Installation", "Number of tracked days",
                                               "Gender: Male","Age","Education  (3 levels)",
                                               "Nationality: German", "Marital status", "Household size",
                                               "Privacy attitudes", "Computer know-how"))


table_descriptives1$N <- colSums(!is.na(gsurf_consentexp_analysis_full[, variables]))
table_descriptives1$Mean <- colMeans(gsurf_consentexp_analysis_full[, variables], na.rm = TRUE)
table_descriptives1$sd <- apply(gsurf_consentexp_analysis_full[, variables], 2, sd, na.rm = TRUE)
table_descriptives1$min <- apply(gsurf_consentexp_analysis_full[, variables], 2, min, na.rm = TRUE)
table_descriptives1$median <- apply(gsurf_consentexp_analysis_full[, variables], 2, median, na.rm = TRUE)
table_descriptives1$max <- apply(gsurf_consentexp_analysis_full[, variables], 2, max, na.rm = TRUE)

stargazer(table_descriptives1, type = "text", 
          style = "asr", align = F, summary=F, rownames = F,
          digits = 3,
          table.placement = "H", no.space = T, single.row = T, header = F, 
          dep.var.caption = NULL)
```



## For table A3 - Meta only
```{r}
table_descriptives2 <- data.frame(Variable = c("Consent", "Installation", "Number of tracked days", "Gender: Male","Age","Education  (3 levels)",
                                               "Nationality: German", "Marital status", "Household size",
                                               "Privacy attitudes", "Computer know-how"))

gsurf_consentexp_analysis_Meta <- gsurf_consentexp_analysis_full %>% filter(rec_arm == "Meta")

table_descriptives2$N <- colSums(!is.na(gsurf_consentexp_analysis_Meta[, variables]))
table_descriptives2$Mean <- colMeans(gsurf_consentexp_analysis_Meta[, variables], na.rm = TRUE)
table_descriptives2$sd <- apply(gsurf_consentexp_analysis_Meta[, variables], 2, sd, na.rm = TRUE)
table_descriptives2$min <- apply(gsurf_consentexp_analysis_Meta[, variables], 2, min, na.rm = TRUE)
table_descriptives2$median <- apply(gsurf_consentexp_analysis_Meta[, variables], 2, median, na.rm = TRUE)
table_descriptives2$max <- apply(gsurf_consentexp_analysis_Meta[, variables], 2, max, na.rm = TRUE)

stargazer(table_descriptives2, type = "text", 
          style = "asr", align = F, summary=F, rownames = F,
          digits = 3,
          table.placement = "H", no.space = T, single.row = T, header = F, 
          dep.var.caption = NULL)
```

## For table A4 - ALLBUS only
```{r}
table_descriptives3 <- data.frame(Variable = c("Consent", "Installation", "Number of tracked days", "Gender: Male","Age","Education  (3 levels)",
                                               "Nationality: German", "Marital status", "Household size",
                                               "Privacy attitudes", "Computer know-how"))

gsurf_consentexp_analysis_ALLBUS <- gsurf_consentexp_analysis_full %>% filter(rec_arm == "ALLBUS")

table_descriptives3$N <- colSums(!is.na(gsurf_consentexp_analysis_ALLBUS[, variables]))
table_descriptives3$Mean <- colMeans(gsurf_consentexp_analysis_ALLBUS[, variables], na.rm = TRUE)
table_descriptives3$sd <- apply(gsurf_consentexp_analysis_ALLBUS[, variables], 2, sd, na.rm = TRUE)
table_descriptives3$min <- apply(gsurf_consentexp_analysis_ALLBUS[, variables], 2, min, na.rm = TRUE)
table_descriptives3$median <- apply(gsurf_consentexp_analysis_ALLBUS[, variables], 2, median, na.rm = TRUE)
table_descriptives3$max <- apply(gsurf_consentexp_analysis_ALLBUS[, variables], 2, max, na.rm = TRUE)

stargazer(table_descriptives3, type = "text", 
          style = "asr", align = F, summary=F, rownames = F,
          digits = 3,
          table.placement = "H", no.space = T, single.row = T, header = F, 
          dep.var.caption = NULL)
```

## For table A5 - Consenters only
```{r}
table_descriptives4 <- data.frame(Variable = c("Consent", "Installation", "Number of tracked days", "Gender: Male","Age","Education  (3 levels)",
                                               "Nationality: German", "Marital status", "Household size",
                                               "Privacy attitudes", "Computer know-how"))

gsurf_consentexp_analysis_consent <- gsurf_consentexp_analysis_full %>% filter(consent == 1)

table_descriptives4$N <- colSums(!is.na(gsurf_consentexp_analysis_consent[, variables]))
table_descriptives4$Mean <- colMeans(gsurf_consentexp_analysis_consent[, variables], na.rm = TRUE)
table_descriptives4$sd <- apply(gsurf_consentexp_analysis_consent[, variables], 2, sd, na.rm = TRUE)
table_descriptives4$min <- apply(gsurf_consentexp_analysis_consent[, variables], 2, min, na.rm = TRUE)
table_descriptives4$median <- apply(gsurf_consentexp_analysis_consent[, variables], 2, median, na.rm = TRUE)
table_descriptives4$max <- apply(gsurf_consentexp_analysis_consent[, variables], 2, max, na.rm = TRUE)

stargazer(table_descriptives4, type = "text", 
          style = "asr", align = F, summary=F, rownames = F, 
          digits = 3,
          table.placement = "H", no.space = T, single.row = T, header = F, 
          dep.var.caption = NULL)
```

## For table A6 - Installers only
```{r}
table_descriptives5 <- data.frame(Variable = c("Consent", "Installation", "Number of tracked days", "Gender: Male","Age","Education  (3 levels)",
                                               "Nationality: German", "Marital status", "Household size",
                                               "Privacy attitudes", "Computer know-how"))

gsurf_consentexp_analysis_installation <- gsurf_consentexp_analysis_full %>% filter(installation == 1)

table_descriptives5$N <- colSums(!is.na(gsurf_consentexp_analysis_installation[, variables]))
table_descriptives5$Mean <- colMeans(gsurf_consentexp_analysis_installation[, variables], na.rm = TRUE)
table_descriptives5$sd <- apply(gsurf_consentexp_analysis_installation[, variables], 2, sd, na.rm = TRUE)
table_descriptives5$min <- apply(gsurf_consentexp_analysis_installation[, variables], 2, min, na.rm = TRUE)
table_descriptives5$median <- apply(gsurf_consentexp_analysis_installation[, variables], 2, median, na.rm = TRUE)
table_descriptives5$max <- apply(gsurf_consentexp_analysis_installation[, variables], 2, max, na.rm = TRUE)

stargazer(table_descriptives5, type = "text", 
          style = "asr", align = F, summary=F, rownames = F, 
          digits = 3,
          table.placement = "H", no.space = T, single.row = T, header = F, 
          dep.var.caption = NULL)
```

## For table A8
```{r}
randomization_test_data <- gsurf_consentexp_analysis_full %>% 
  mutate(group = if_else(uncon_inc == 0 & con_inc == 1,1,
                         if_else(uncon_inc == 0 & con_inc == 2,2,
                                 if_else(uncon_inc == 0 & con_inc == 3,3,
                                         if_else(uncon_inc == 0 & con_inc == 4,4,
                                                 if_else(uncon_inc == 1 & con_inc == 1,5,
                                                         if_else(uncon_inc == 1 & con_inc == 2,6,
                                                                 if_else(uncon_inc == 1 & con_inc == 3,7,
                                                                         if_else(uncon_inc == 1 & con_inc == 4,8,0)))))))))


sexX <- chisq.test(table(randomization_test_data$group,randomization_test_data$sex))
germanX <- chisq.test(table(randomization_test_data$group,randomization_test_data$german))
marriedX <- chisq.test(table(randomization_test_data$group,randomization_test_data$married))
educX <- chisq.test(table(randomization_test_data$group,randomization_test_data$educ))

df_chisq <- data.frame(demographic = c("Gender: Male","Nationality: German","Marital status","Education (3 levels)"),
                       test_statistic = round(c(sexX$statistic, germanX$statistic, marriedX$statistic, educX$statistic), digits = 3),
                       pvalue = round(c(sexX$p.value, germanX$p.value, marriedX$p.value, educX$p.value), digits = 3))

#ANOVA
ageA <- oneway.test(age~  group,
  data = randomization_test_data, var.equal = T)

hhsizA <- oneway.test(hhsiz~  group,
  data = randomization_test_data, var.equal = T)

paA <- oneway.test(pa_index~  group,
  data = randomization_test_data, var.equal = T)

khA <-oneway.test(kh_index~  group,
  data = randomization_test_data, var.equal = T)

df_anova <- data.frame(demographic = c("Age","Household size", "Privacy attitudes", "Computer know-how"),
                       test_statistic = round(c(ageA$statistic, hhsizA$statistic, paA$statistic, khA$statistic), digits = 3),
                       pvalue = round(c(ageA$p.value, hhsizA$p.value, paA$p.value, khA$p.value), digits = 3))

df_randomization <- rbind(df_chisq,df_anova)
stargazer(df_randomization, type = "text", 
          style = "asr", align = F, summary=F, rownames = F, 
          digits = 3,
          table.placement = "H", no.space = T, single.row = T, header = F, 
          dep.var.caption = NULL)
```

# for table A9 (those who gave consent and those who installed)
```{r}
# Chi-square consent
randomization_test_data_consent <- gsurf_consentexp_analysis_full %>% 
  filter(consent == 1) %>% 
  mutate(group = if_else(uncon_inc == 0 & con_inc == 1,1,
                         if_else(uncon_inc == 0 & con_inc == 2,2,
                                 if_else(uncon_inc == 0 & con_inc == 3,3,
                                         if_else(uncon_inc == 0 & con_inc == 4,4,
                                                 if_else(uncon_inc == 1 & con_inc == 1,5,
                                                         if_else(uncon_inc == 1 & con_inc == 2,6,
                                                                 if_else(uncon_inc == 1 & con_inc == 3,7,
                                                                         if_else(uncon_inc == 1 & con_inc == 4,8,0)))))))))


sexX <- chisq.test(table(randomization_test_data_consent$group,randomization_test_data_consent$sex))
germanX <- chisq.test(table(randomization_test_data_consent$group,randomization_test_data_consent$german))
marriedX <- chisq.test(table(randomization_test_data_consent$group,randomization_test_data_consent$married))
educX <- chisq.test(table(randomization_test_data_consent$group,randomization_test_data_consent$educ))

df_chisq_consent <- data.frame(demographic = c("Sex","German nationality","Maritial status: Married","Education (3 levels)"),
                            test_statistic = round(c(sexX$statistic, germanX$statistic, marriedX$statistic, educX$statistic), digits = 3),
                            pvalue = round(c(sexX$p.value, germanX$p.value, marriedX$p.value, educX$p.value), digits = 3))

#Chi-square installation
randomization_test_data_installation <- gsurf_consentexp_analysis_full %>% 
  filter(installation == 1) %>% 
  mutate(group = if_else(uncon_inc == 0 & con_inc == 1,1,
                         if_else(uncon_inc == 0 & con_inc == 2,2,
                                 if_else(uncon_inc == 0 & con_inc == 3,3,
                                         if_else(uncon_inc == 0 & con_inc == 4,4,
                                                 if_else(uncon_inc == 1 & con_inc == 1,5,
                                                         if_else(uncon_inc == 1 & con_inc == 2,6,
                                                                 if_else(uncon_inc == 1 & con_inc == 3,7,
                                                                         if_else(uncon_inc == 1 & con_inc == 4,8,0)))))))))


sexX <- chisq.test(table(randomization_test_data_installation$group,randomization_test_data_installation$sex))
germanX <- chisq.test(table(randomization_test_data_installation$group,randomization_test_data_installation$german))
marriedX <- chisq.test(table(randomization_test_data_installation$group,randomization_test_data_installation$married))
educX <- chisq.test(table(randomization_test_data_installation$group,randomization_test_data_installation$educ))

df_chisq_installation <- data.frame(demographic = c("Sex","German nationality","Maritial status: Married","Education (3 levels)"),
                              test_statistic = round(c(sexX$statistic, germanX$statistic, marriedX$statistic, educX$statistic), digits = 3),
                              pvalue = round(c(sexX$p.value, germanX$p.value, marriedX$p.value, educX$p.value), digits = 3))

#joined output table
df_chisq_coninst <- left_join(df_chisq_consent, df_chisq_installation, by = "demographic", suffix = c(".consent", ".installation"))


#ANOVA

#consent
ageA_consent <- oneway.test(age~  group,
                    data = randomization_test_data_consent, var.equal = T)

hhsizA <- oneway.test(hhsiz~  group,
                      data = randomization_test_data_consent, var.equal = T)

paA <- oneway.test(pa_index~  group,
                   data = randomization_test_data_consent, var.equal = T)

khA <-oneway.test(kh_index~  group,
                  data = randomization_test_data_consent, var.equal = T)

df_anova_consent <- data.frame(demographic = c("Age","Household size","Privacy attitudes","Computer know-how"),
                            test_statistic = round(c(ageA_consent$statistic, hhsizA$statistic, paA$statistic, khA$statistic), digits = 3),
                            pvalue = round(c(ageA_consent$p.value, hhsizA$p.value, paA$p.value, khA$p.value), digits = 3))
#installation
ageA_installation <- oneway.test(age~  group,
                    data = randomization_test_data_installation, var.equal = T)

hhsizA <- oneway.test(hhsiz~  group,
                      data = randomization_test_data_installation, var.equal = T)

paA <- oneway.test(pa_index~  group,
                   data = randomization_test_data_installation, var.equal = T)

khA <-oneway.test(kh_index~  group,
                  data = randomization_test_data_installation, var.equal = T)

df_anova_installation <- data.frame(demographic = c("Age","Household size","Privacy attitudes","Computer know-how"),
                              test_statistic = round(c(ageA_installation$statistic, hhsizA$statistic, paA$statistic, khA$statistic), digits = 3),
                              pvalue = round(c(ageA_installation$p.value, hhsizA$p.value, paA$p.value, khA$p.value), digits = 3))

#joined output table
df_anova_coninst <- left_join(df_anova_consent, df_anova_installation, by = "demographic", suffix = c(".consent", ".installation"))


df_randomization_coninst <- rbind(df_chisq_coninst,df_anova_coninst)
stargazer(df_randomization_coninst, type = "text", 
          style = "asr", align = F, summary=F, rownames = F, 
          digits = 3,
          table.placement = "H", no.space = T, single.row = T, header = F, 
          dep.var.caption = NULL)
```

## Table A10: Mean Age and Standard Deviations Across All Eight Incentive Conditions

```{r}
randomization_age <- gsurf_consentexp_analysis_full %>% 
  group_by(uncon_inc, con_inc) %>% 
  summarise(mean_age_full_sample = as.character(round(mean(age, na.rm = T),3)),
            sd_age_full_sample = as.character(round(sd(age, na.rm = T),3))) %>% 
  left_join(gsurf_consentexp_analysis_full %>% 
              filter(consent == 1) %>% 
              group_by(uncon_inc, con_inc) %>% 
              summarise(mean_age_consent_given = as.character(round(mean(age, na.rm = T),3)),
                        sd_age_consent_given = as.character(round(sd(age, na.rm = T),3))), by = c("uncon_inc","con_inc")) %>% 
  left_join(gsurf_consentexp_analysis_full %>% 
              filter(installation == 1) %>% 
              group_by(uncon_inc, con_inc) %>% 
              summarise(mean_age_installed = as.character(round(mean(age, na.rm = T),3)),
                        sd_age_installed = as.character(round(sd(age, na.rm = T),3))), by = c("uncon_inc","con_inc")) %>% 
  mutate(uncon_inc = case_match(uncon_inc,
                                0 ~ "none",
                                1 ~ "5€"),
         con_inc = case_match(con_inc,
                              1 ~ "none",
                              2 ~ "10€",
                              3 ~ "25€",
                              4 ~ "40€"))
randomization_age_tab <- randomization_age %>% 
  mutate(mean_age_full_sample = paste0(mean_age_full_sample, " (", sd_age_full_sample, ")"),
         mean_age_consent_given = paste0(mean_age_consent_given, " (", sd_age_consent_given, ")"),
         mean_age_installed = paste0(mean_age_installed, " (", sd_age_installed, ")")) %>% 
  select(uncon_inc,con_inc,mean_age_full_sample,mean_age_consent_given,mean_age_installed)


stargazer(randomization_age_tab, type = "text", style = "asr", align = F, summary=F, rownames = F, covariate.labels = c("Uncon. Inc.", "Con. Inc.", "Mean (SD) Full Sample", "Mean (SD) Consent Given", "Mean (SD) Installed"), digits = 1, table.placement = "H", no.space = T, single.row = T, header = F, 
          dep.var.caption = NULL)
```
## For text in Subsection A3 Randomization checks

```{r}
cat("SD of the condition means in age for consenters:")
sd(randomization_age$mean_age_consent_given)
cat("SD of the condition means in age for installers:")
sd(randomization_age$mean_age_installed)
cat("\n")

cat("Overall SD of age among consenters:")
sd(gsurf_consentexp_analysis_consent$age)
cat("Overall SD of age among installers:")
sd(gsurf_consentexp_analysis_installation$age)
cat("\n")

cat("T-Test for age grouped by unconditional/ no unconditional incentive among consenters:")
t.test(randomization_test_data_consent$age ~ randomization_test_data_consent$uncon_inc)
cat("T-Test for age grouped by conditional/ no conditional incentive among consenters:")
t.test(randomization_test_data_consent$age ~ randomization_test_data_consent$con_inc_total)
cat("\n")

cat("T-Test for age grouped by unconditional/ no unconditional incentive among installers:")
t.test(randomization_test_data_installation$age ~ randomization_test_data_installation$uncon_inc)
cat("T-Test for age grouped by conditional/ no conditional incentive among installers:")
t.test(randomization_test_data_installation$age ~ randomization_test_data_installation$con_inc_total)
```

### footnote 5

```{r}
model_group_age_consent <- lmer(age ~ 1 + (1 | group), data = randomization_test_data_consent)
icc(model_group_age_consent)

model_group_age_installation <- lmer(age ~ 1 + (1 | group), data = randomization_test_data_installation)
icc(model_group_age_installation)
```


## For table A11 - Unadjusted Proportion and Means by Incentive Condition
unadjusted results

```{r}
table_unadjusted <- data.frame(Condition = c("no unconditional","5€ unconditional","no conditional","10€ conditional","25€ conditional","40€ conditional"))

table_unadjusted <- table_unadjusted %>% 
  mutate(Consent = c(((gsurf_consentexp_analysis_full %>%  group_by(uncon_inc) %>% summarize(consent = mean(consent)))$consent), ((gsurf_consentexp_analysis_full %>%  group_by(con_inc) %>% summarize(consent = mean(consent)))$consent)),
         Installation = c(((gsurf_consentexp_analysis_full %>%  group_by(uncon_inc) %>% summarize(installation = mean(installation, na.rm = T)))$installation), ((gsurf_consentexp_analysis_full %>%  group_by(con_inc) %>% summarize(installation = mean(installation, na.rm = T)))$installation)),
         Tracked_Days = c(((gsurf_consentexp_analysis_full %>%  group_by(uncon_inc) %>% summarize(active_days_threshold = mean(active_days_threshold, na.rm = T)))$active_days_threshold), ((gsurf_consentexp_analysis_full %>%  group_by(con_inc) %>% summarize(active_days_threshold = mean(active_days_threshold, na.rm = T)))$active_days_threshold)))

stargazer(table_unadjusted, type = "text", 
          style = "asr", align = F, summary=F, rownames = F, #covariate.labels = c(),
          digits = 3,
          table.placement = "H", no.space = T, single.row = T, header = F, 
          dep.var.caption = NULL)

```

# Model output

## Table A12: Logistic Regression for Outcome Variable Consent (Full Output)

```{r, warning=FALSE}
glm_consent_pa_kh_demo_coninc_dummy <- glm(consent ~ uncon_inc+con_inc_total+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full, family = "binomial")

glm_consent_pa_kh_demo_dummy <- glm(consent ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full, family = "binomial")

labels <- c("Unconditional 5€","Conditional any amount","Conditional 10€","Conditional 25€","Conditional 40€",
            "Privacy attitudes","Computer know-how","Age","Male","Education: medium (ref.: low)",
            "Education: high (ref.: low)","Married","Household size", "German nationality",
            "Recruitment via Meta (ref.: ALLBUS)","Study number (ref.: study two)", "Intercept")

stargazer(glm_consent_pa_kh_demo_coninc_dummy, glm_consent_pa_kh_demo_dummy, type = "text", style = "default", align = F, summary=T, rownames = T, covariate.labels = labels, table.placement = "H", no.space = T, single.row = T, header = F, star.cutoffs = NA,
          dep.var.caption = NULL, report = "vcstp", dep.var.labels = "Consent")

```

## Logit marginal effects consent

```{r, warning=FALSE}
## Unconditional incentive

# Get estimated marginal means (EMMs)
emm_consent_uncon_inc <- emmeans(glm_consent_pa_kh_demo_dummy, ~ uncon_inc, type = "response", weights = "proportional")
summary(emm_consent_uncon_inc, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_consent_uncon_inc, method = "pairwise", adjust = "tukey")
summary(contrast_results)

## Conditional incentive dummy

# Get estimated marginal means (EMMs)
emm_consent_con_inc_dummy <- emmeans(glm_consent_pa_kh_demo_coninc_dummy, ~ con_inc_total, type = "response", weights = "proportional")
summary(emm_consent_con_inc_dummy, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_consent_con_inc_dummy, method = "pairwise", adjust = "tukey")
summary(contrast_results)

## Conditional incentive

# Get estimated marginal means (EMMs)
emm_consent_con_inc <- emmeans(glm_consent_pa_kh_demo_dummy, ~ con_inc, type = "response", weights = "proportional")
summary(emm_consent_con_inc, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_consent_con_inc, method = "pairwise", adjust = "tukey")
summary(contrast_results)

## Recruitment method
# Get estimated marginal means (EMMs)
emm_consent_rec_arm <- emmeans(glm_consent_pa_kh_demo_dummy, ~ rec_arm, type = "response", weights = "proportional")
summary(emm_consent_rec_arm, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_consent_rec_arm, method = "pairwise", adjust = "tukey")
summary(contrast_results)
```


## Table A13: Logistic Regression for Outcome Variable Installation (Full Output)

```{r, warning=FALSE}
glm_installation_pa_kh_demo_coninc_dummy_consent <- glm(installation ~ uncon_inc+con_inc_total+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full_consent, family = "binomial")

glm_installation_pa_kh_demo_dummy_consent <- glm(installation ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full_consent, family = "binomial")

labels <- c("Unconditional 5€","Conditional any amount","Conditional 10€","Conditional 25€","Conditional 40€",
            "Privacy attitudes","Computer know-how","Age","Male","Education: medium (ref.: low)",
            "Education: high (ref.: low)","Married","Household size", "German nationality",
            "Recruitment via Meta (ref.: ALLBUS)","Study number (ref.: study two)", "Intercept")

stargazer(glm_installation_pa_kh_demo_coninc_dummy_consent, glm_installation_pa_kh_demo_dummy_consent, type = "text", style = "default", align = F, summary=T, rownames = T, covariate.labels = labels, table.placement = "H", no.space = T, single.row = T, header = F, star.cutoffs = NA,
          dep.var.caption = NULL, dep.var.labels = "Installation", report = "vcstp")

```


## Logit marginal effects installation
```{r, warning=FALSE}
## Unconditional incentive

# Get estimated marginal means (EMMs)
emm_installation_uncon_inc <- emmeans(glm_installation_pa_kh_demo_dummy_consent, ~ uncon_inc, type = "response", weights = "proportional")
summary(emm_installation_uncon_inc, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_installation_uncon_inc, method = "pairwise", adjust = "tukey")
summary(contrast_results)

## Conditional incentive dummy

# Get estimated marginal means (EMMs)
emm_installation_con_inc_dummy <- emmeans(glm_installation_pa_kh_demo_coninc_dummy_consent, ~ con_inc_total, type = "response", weights = "proportional")
summary(emm_installation_con_inc_dummy, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_installation_con_inc_dummy, method = "pairwise", adjust = "tukey")
summary(contrast_results)

## Conditional incentive

# Get estimated marginal means (EMMs)
emm_installation_con_inc <- emmeans(glm_installation_pa_kh_demo_dummy_consent, ~ con_inc, type = "response", weights = "proportional")
summary(emm_installation_con_inc, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_installation_con_inc, method = "pairwise", adjust = "tukey")
summary(contrast_results)

## Recruitment method
# Get estimated marginal means (EMMs)
emm_installation_rec_arm <- emmeans(glm_installation_pa_kh_demo_dummy_consent, ~ rec_arm, type = "response", weights = "proportional")
summary(emm_installation_rec_arm, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_installation_rec_arm, method = "pairwise", adjust = "tukey")
summary(contrast_results)
```

## Table A14: Negative Binomial Regression for Outcome Number of Tracked Days (Full Output)

```{r}
lm_active_full_transfer_coninc_dummy <- glmmTMB(active_days_threshold  ~ uncon_inc+con_inc_total+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, family=nbinom2, data = gsurf_consentexp_analysis_full_active)

lm_active_full_transfer_demo <- glmmTMB(active_days_threshold  ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, family=nbinom2, data = gsurf_consentexp_analysis_full_active)
nobs(lm_active_full_transfer_demo)
```


testing for overdispersion
```{r}
check_overdispersion(lm_active_full_transfer_coninc_dummy, alternative = c("two.sided"))
check_overdispersion(lm_active_full_transfer_demo, alternative = c("two.sided"))
```

```{r, include = FALSE}
labels <- c("Unconditional 5€","","","Conditional any amount","","","Conditional 10€","","","Conditional 25€","","","Conditional 40€","","",
            "Privacy attitudes","","","Computer know-how","","","Age","","","Male","","","Education: medium (ref.: low)","","",
            "Education: high (ref.: low)","","","Married","","","Household size","","", "German nationality","","",
            "Recruitment via Meta (ref.: ALLBUS)","","","Study number (ref.: study two)","","","Intercept","","")

model1a <- data.frame(Estimate = unlist(lapply(coef(summary(lm_active_full_transfer_coninc_dummy)),function(x) x[,"Estimate"])),
                      stderr = paste0("(",as.character(round(unlist(lapply(coef(summary(lm_active_full_transfer_coninc_dummy)),function(x) x[,"Std. Error"])),3)),")"),
                      zvalue = unlist(lapply(coef(summary(lm_active_full_transfer_coninc_dummy)),function(x) x[,"z value"])),
                      pvalue = unlist(lapply(coef(summary(lm_active_full_transfer_coninc_dummy)),function(x) x[,"Pr(>|z|)"]))) %>%
  mutate(Effect = rownames(.),
         p_value = if_else(pvalue <0.001, "p < 0.001",
                         paste0("p = ",paste0(sprintf("%.3f",round(pvalue, digits = 3))))),
         z_value = paste0("z = ",paste0(sprintf("%.3f",round(zvalue, digits = 3)))),
         Effectestimate = paste0(sprintf("%.3f",round(Estimate, digits = 3))," ",stderr)) %>% 
  select(Effect, Effectestimate, z_value, p_value) %>% 
  pivot_longer(-c(Effect), names_to = "statistic", values_to = "value")

model1b <- data.frame(Estimate = unlist(lapply(coef(summary(lm_active_full_transfer_demo)),function(x) x[,"Estimate"])),
                      stderr = paste0("(",as.character(round(unlist(lapply(coef(summary(lm_active_full_transfer_demo)),function(x) x[,"Std. Error"])),3)),")"),
                      zvalue = unlist(lapply(coef(summary(lm_active_full_transfer_demo)),function(x) x[,"z value"])),
                      pvalue = unlist(lapply(coef(summary(lm_active_full_transfer_demo)),function(x) x[,"Pr(>|z|)"]))) %>% 
  mutate(Effect = rownames(.),
         p_value = if_else(pvalue <0.001, "p < 0.001",
                         paste0("p = ",paste0(sprintf("%.3f",round(pvalue, digits = 3))))),
         z_value = paste0("z = ",paste0(sprintf("%.3f",round(zvalue, digits = 3)))),
         Effectestimate = paste0(sprintf("%.3f",round(Estimate, digits = 3))," ",stderr)) %>% 
  select(Effect, Effectestimate, z_value, p_value) %>% 
  pivot_longer(-c(Effect), names_to = "statistic", values_to = "value")




models <- full_join(model1a,model1b, by = c("Effect", "statistic")) %>% 
  mutate(n = c(49:51,1:6,16:48,7:15)) %>% 
  arrange(n) %>% 
  mutate(Effects = labels) %>% 
  select(-Effect, -n, -statistic) %>% 
  select(Effects, model1 = value.x, model2 =value.y) %>% 
  rbind(c("N",as.character(nobs(lm_active_full_transfer_coninc_dummy)),as.character(nobs(lm_active_full_transfer_demo))),
        c("Log Likelihood", as.character(round(summary(lm_active_full_transfer_coninc_dummy)$AIC[3],digits = 3)),as.character(round(summary(lm_active_full_transfer_demo)$AIC[3],digits = 3))),
        c("AIC", as.character(round(summary(lm_active_full_transfer_coninc_dummy)$AIC[1],digits = 3)),as.character(round(summary(lm_active_full_transfer_demo)$AIC[1],digits = 3))))
```

```{r}
stargazer(models, type = "text", style = "asr", align = F, summary=F, rownames = F, covariate.labels = c("","(1)","(2)"), table.placement = "H", no.space = T, single.row = T, header = F, dep.var.caption = NULL, title = "Tracked Days")
```

## Average Marginal Effects Tracked Days

```{r}
## Unconditional incentive

# Get estimated marginal means (EMMs)
emm_activdays_uncon_inc <- emmeans(lm_active_full_transfer_demo, ~ uncon_inc, type = "response", weights = "proportional")
summary(emm_activdays_uncon_inc, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_activdays_uncon_inc, method = "pairwise", adjust = "tukey") 
summary(contrast_results)

## Conditional incentive dummy

# Get estimated marginal means (EMMs)
emm_activdays_con_inc_dummy <- emmeans(lm_active_full_transfer_coninc_dummy, ~ con_inc_total, type = "response", weights = "proportional")
summary(emm_activdays_con_inc_dummy, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_activdays_con_inc_dummy, method = "pairwise", adjust = "tukey") 
summary(contrast_results)

## Conditional incentive

# Get estimated marginal means (EMMs)
emm_activdays_con_inc <- emmeans(lm_active_full_transfer_demo, ~ con_inc, type = "response", weights = "proportional")
summary(emm_activdays_con_inc, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_activdays_con_inc, method = "pairwise", adjust = "tukey") 
summary(contrast_results)

## Recruitment method
# Get estimated marginal means (EMMs)
emm_activdays_rec_arm <- emmeans(lm_active_full_transfer_demo, ~ rec_arm, type = "response", weights = "proportional")
summary(emm_activdays_rec_arm, level = 0.834) #for CIs at 83.4% intervalls
# Pairwise comparisons with p-values
contrast_results <- contrast(emm_activdays_rec_arm, method = "pairwise", adjust = "tukey") 
summary(contrast_results)
```

## Table A15: Logistic Regression for Outcome Variable Consent Considering Interactions Between Incentives and Recruitment Methods (Full Output)

```{r}
glm_consent_recarm_unconinc_interact <- glm(consent ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num+rec_arm*uncon_inc, data = gsurf_consentexp_analysis_full, family = "binomial")

glm_consent_recarm_coninc_interact <- glm(consent ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num+rec_arm*as.factor(con_inc), data = gsurf_consentexp_analysis_full, family = "binomial")
```

```{r}
labels_interaction <- c("Unconditional 5€","Conditional 10€","Conditional 25€","Conditional 40€",
            "Privacy attitudes","Computer know-how","Age","Sex","Education: medium (ref.: low)",
            "Education: high (ref.: low)","Married","Household size", "German nationality",
            "Recruitment via Meta (ref.: ALLBUS)","Study number (ref.: study two)","Unconditional 5€ X recruitment via Meta",
            "Conditional 10€ X recruitment via Meta","Conditional 25€ X recruitment via Meta","Conditional 40€ X recruitment via Meta", "Intercept")

stargazer(glm_consent_recarm_unconinc_interact, glm_consent_recarm_coninc_interact, type = "text", style = "default", align = F, summary=T, rownames = T, covariate.labels = labels_interaction, table.placement = "H", no.space = T, single.row = T, header = F, star.cutoffs = NA, dep.var.caption = NULL, report = "vcstp", dep.var.labels = "Consent")
```

## Table A16: Logistic Regression for Outcome Variable Installation Considering Interactions Between Incentives and Recruitment Methods (Full Output)

```{r}
glm_installation_recarm_unconinc_interact <- glm(installation ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num+rec_arm*uncon_inc, data = gsurf_consentexp_analysis_full, family = "binomial")

glm_installation_recarm_coninc_interact <- glm(installation ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num+rec_arm*as.factor(con_inc), data = gsurf_consentexp_analysis_full, family = "binomial")
```

```{r}
labels_interaction <- c("Unconditional 5€","Conditional 10€","Conditional 25€","Conditional 40€",
            "Privacy attitudes","Computer know-how","Age","Sex","Education: medium (ref.: low)",
            "Education: high (ref.: low)","Married","Household size", "German nationality",
            "Recruitment via Meta (ref.: ALLBUS)","Study number (ref.: study two)","Unconditional 5€ X recruitment via Meta",
            "Conditional 10€ X recruitment via Meta","Conditional 25€ X recruitment via Meta","Conditional 40€ X recruitment via Meta", "Intercept")

stargazer(glm_installation_recarm_unconinc_interact, glm_installation_recarm_coninc_interact, type = "text", style = "default", align = F, summary=T, rownames = T, covariate.labels = labels_interaction, table.placement = "H", no.space = T, single.row = T, header = F, star.cutoffs = NA, dep.var.caption = NULL, report = "vcstp", dep.var.labels = "Installation")
```

## Table A17: Negative Binomial Regression for Outcome Variable Number of Tracked Days Considering Interactions Between Incentives and Recruitment Methods (Full Output)

```{r}
lm_active_full_transfer_demo_unconinc_interact <- glmmTMB(active_days_threshold  ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num+rec_arm*uncon_inc, family=nbinom2, data = gsurf_consentexp_analysis_full)

lm_active_full_transfer_demo_coninc_interact <- glmmTMB(active_days_threshold  ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num+rec_arm*as.factor(con_inc), family=nbinom2, data = gsurf_consentexp_analysis_full)
```


```{r, include=FALSE}
labels_interaction <- c("Unconditional 5€","","","Conditional 10€","","","Conditional 25€","","","Conditional 40€","","",
            "Privacy attitudes","","","Computer know-how","","","Age","","","Male","","","Education: medium (ref.: low)","","",
            "Education: high (ref.: low)","","","Married","","","Household size","","", "German nationality","","",
            "Recruitment via Meta (ref.: ALLBUS)","","","Study number (ref.: study two)","","","Unconditional 5€ X recruitment via Meta","","",
                        "Conditional 10€ X recruitment via Meta","","","conditional 25€ X recruitment via Meta","","",
            "conditional 40€ X recruitment via Meta","","","Intercept","","")


model1a <- data.frame(Estimate = unlist(lapply(coef(summary(lm_active_full_transfer_demo_unconinc_interact)),function(x) x[,"Estimate"])),
                      stderr = paste0("(",as.character(round(unlist(lapply(coef(summary(lm_active_full_transfer_demo_unconinc_interact)),function(x) x[,"Std. Error"])),3)),")"),
                      zvalue = unlist(lapply(coef(summary(lm_active_full_transfer_demo_unconinc_interact)),function(x) x[,"z value"])),
                      pvalue = unlist(lapply(coef(summary(lm_active_full_transfer_demo_unconinc_interact)),function(x) x[,"Pr(>|z|)"]))) %>%
  mutate(Effect = rownames(.),
         p_value = if_else(pvalue <0.001, "p < 0.001",
                           paste0("p = ",paste0(sprintf("%.3f",round(pvalue, digits = 3))))),
         z_value = paste0("z = ",paste0(sprintf("%.3f",round(zvalue, digits = 3)))),
         Effectestimate = paste0(sprintf("%.3f",round(Estimate, digits = 3))," ",stderr)) %>% 
  select(Effect, Effectestimate, z_value, p_value) %>% 
  pivot_longer(-c(Effect), names_to = "statistic", values_to = "value")

model1b <- data.frame(Estimate = unlist(lapply(coef(summary(lm_active_full_transfer_demo_coninc_interact)),function(x) x[,"Estimate"])),
                      stderr = paste0("(",as.character(round(unlist(lapply(coef(summary(lm_active_full_transfer_demo_coninc_interact)),function(x) x[,"Std. Error"])),3)),")"),
                      zvalue = unlist(lapply(coef(summary(lm_active_full_transfer_demo_coninc_interact)),function(x) x[,"z value"])),
                      pvalue = unlist(lapply(coef(summary(lm_active_full_transfer_demo_coninc_interact)),function(x) x[,"Pr(>|z|)"]))) %>% 
   mutate(Effect = rownames(.),
         p_value = if_else(pvalue <0.001, "p < 0.001",
                           paste0("p = ",paste0(sprintf("%.3f",round(pvalue, digits = 3))))),
         z_value = paste0("z = ",paste0(sprintf("%.3f",round(zvalue, digits = 3)))),
         Effectestimate = paste0(sprintf("%.3f",round(Estimate, digits = 3))," ",stderr)) %>% 
  select(Effect, Effectestimate, z_value, p_value) %>% 
  pivot_longer(-c(Effect), names_to = "statistic", values_to = "value")



models <- full_join(model1a,model1b, by = c("Effect", "statistic")) %>% 
  mutate(n = c(58:60,1:57)) %>% 
  arrange(n) %>% 
  mutate(Effects = labels_interaction) %>% 
  select(-Effect, -n, -statistic) %>% 
  select(Effects, model1 = value.x, model2 =value.y) %>% 
  rbind(c("N",as.character(nobs(lm_active_full_transfer_demo_unconinc_interact)),as.character(nobs(lm_active_full_transfer_demo_coninc_interact))),
        c("Log Likelihood", as.character(round(summary(lm_active_full_transfer_demo_unconinc_interact)$AIC[3],digits = 3)),as.character(round(summary(lm_active_full_transfer_demo_coninc_interact)$AIC[3],digits = 3))),
        c("AIC", as.character(round(summary(lm_active_full_transfer_demo_unconinc_interact)$AIC[1],digits = 3)),as.character(round(summary(lm_active_full_transfer_demo_coninc_interact)$AIC[1],digits = 3))))


```

```{r}
stargazer(models, type = "text", style = "asr", align = F, summary=F, rownames = F, covariate.labels = c("","(1)","(2)"), table.placement = "H", no.space = T, single.row = T,
          header = F, dep.var.caption = NULL, dep.var.labels = "Tracked days", title = "Tracked Days")
```

## Table 18: Linear Probability Models for the Outcome Variable Consent

```{r}
lm_consent_pa_kh_demo_dummy <- lm(consent ~ uncon_inc+con_inc_total+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full)

lm_consent_pa_kh_demo <- lm(consent ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full)
```

```{r}
labels <- c("Unconditional 5€","Conditional any amount","Conditional 10€","Conditional 25€","Conditional 40€",
            "Privacy attitudes","Computer know-how","Age","Male","Education: medium (ref.: low)",
            "Education: high (ref.: low)","Married","Household size", "German nationality",
            "Recruitment via Meta (ref.: ALLBUS)","Study number (ref.: study two)", "Intercept")

stargazer(lm_consent_pa_kh_demo_dummy, lm_consent_pa_kh_demo,  type = "text", style = "default", align = F, summary=T, rownames = T, covariate.labels = labels, table.placement = "H", no.space = T, single.row = T, header = F, star.cutoffs = NA, dep.var.caption = NULL, report = "vcstp", dep.var.labels = "Consent")
```

## Table A19: Linear Probability Models for the Outcome Variable Installation

```{r}
lm_installation_pa_kh_demo_dummy <- lm(installation ~ uncon_inc+ con_inc_total+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full_consent)

lm_installation_pa_kh_demo <- lm(installation ~ uncon_inc+ as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full_consent)
```

```{r}
labels <- c("Unconditional 5€","Conditional any amount","Conditional 10€","Conditional 25€","Conditional 40€",
            "Privacy attitudes","Computer know-how","Age","Male","Education: medium (ref.: low)",
            "Education: high (ref.: low)","Married","Household size", "German nationality",
            "Recruitment via Meta (ref.: ALLBUS)","Study number (ref.: study two)","Intercept")

stargazer(lm_installation_pa_kh_demo_dummy, lm_installation_pa_kh_demo,  type = "text", style = "default", align = F, summary=T, rownames = T, covariate.labels = labels, table.placement = "H", no.space = T, single.row = T, header = F, star.cutoffs = NA, dep.var.caption = NULL, report = "vcstp", dep.var.labels = "Installation")
```


## Table A20: Logistic Regression Models for the Outcome Variable Consent with Bootstrapped Standard Errors

```{r, include=FALSE}
# Define a function for bootstrapping
boot_fn <- function(data, indices, response, predictors) {
  # Take a bootstrap sample
  d <- data[indices, ]
  # Create the formula for the model
  formula <- as.formula(paste(response, "~", paste(predictors, collapse = " + ")))
  # Fit the model on the bootstrap sample
  fit <- glm(formula, data = d, family = binomial)
  # Return the coefficients
  return(coef(fit))
}
```

```{r}
response <- "consent"
set.seed(123)

## conditional incentives binary

predictors_total_con_inc <- c("uncon_inc", "con_inc_total","pa_index","kh_index","age","sex","educ_medium","educ_high","married","hhsiz","german","rec_arm","study_tracking_num")

mymodel_consent_total_con_inc <- glm_consent_pa_kh_demo_coninc_dummy

boot_results_consent_total_con_inc <- boot(data = gsurf_consentexp_analysis_full, statistic = function(data, indices) boot_fn(gsurf_consentexp_analysis_full, indices, response, predictors_total_con_inc), R = 1000)

# Extract bootstrapped standard errors
boot_se <- apply(boot_results_consent_total_con_inc$t, 2, sd)

# Display the original model coefficients with bootstrapped standard errors
original_coef <- coef(mymodel_consent_total_con_inc)
# Calculate z-values and p-values
z_values <- original_coef / boot_se
p_values <- round(2 * pnorm(-abs(z_values)), digits = 3)

# Create a results data frame
results_consent_total_con_inc <- data.frame(
  Estimate = original_coef,
  `Boot SE` = boot_se,
  `z value` = z_values,
  `Pr(>|z|)` = p_values
)

print(results_consent_total_con_inc)

## conditional incentives as factor
predictors_factor_con_inc <- c("uncon_inc", "as.factor(con_inc)","pa_index","kh_index","age","sex","educ_medium","educ_high","married","hhsiz","german","rec_arm","study_tracking_num")

mymodel_consent_factor_con_inc <- glm_consent_pa_kh_demo_dummy


boot_results_consent_factor_con_inc <- boot(data = gsurf_consentexp_analysis_full, statistic = function(data, indices) boot_fn(gsurf_consentexp_analysis_full, indices, response, predictors_factor_con_inc), R = 1000)

# Extract bootstrapped standard errors
boot_se <- apply(boot_results_consent_factor_con_inc$t, 2, sd)

# Display the original model coefficients with bootstrapped standard errors
original_coef <- coef(mymodel_consent_factor_con_inc)
# Calculate z-values and p-values
z_values <- original_coef / boot_se
p_values <- round(2 * pnorm(-abs(z_values)), digits = 3)

# Create a results data frame
results_consent_factor_con_inc <- data.frame(
  Estimate = original_coef,
  `Boot SE` = boot_se,
  `z value` = z_values,
  `Pr(>|z|)` = p_values
)

#print(results_consent_factor_con_inc)
```

```{r, include=FALSE}
labels_boot <- c("Unconditional 5€","","","Conditional any amount","","","Conditional 10€","","","Conditional 25€","","","Conditional 40€","","",
            "Privacy attitudes","","","Computer know-how","","","Age","","","Male","","","Education: medium (ref.: low)","","",
            "Education: high (ref.: low)","","","Married","","","Household size","", "","German nationality","","",
            "Recruitment via Meta (ref.: ALLBUS)","","","Study number (ref.: study two)","","", "Intercept","","")

model1 <- results_consent_total_con_inc %>% 
  mutate(p_value = if_else(`Pr...z..` <0.001, "p < 0.001",
                           paste0("p = ",paste0(sprintf("%.3f",round(`Pr...z..`, digits = 3))))),
         z_value = paste0("z = ",paste0(sprintf("%.3f",round(`z.value`, digits = 3)))),
         stderr = paste0("(",as.character(format(round(`Boot.SE`, digits = 3),nsmall =3)),")"),
         Effectestimate = paste0(sprintf("%.3f",round(Estimate, digits = 3))," ",stderr),
         Effect = rownames(.)) %>% 
  select(Effect, Effectestimate, z_value, p_value) %>%
  pivot_longer(-c(Effect), names_to = "statistic", values_to = "value")



model2 <- results_consent_factor_con_inc %>% 
  mutate(p_value = if_else(`Pr...z..` <0.001, "p < 0.001",
                           paste0("p = ",paste0(sprintf("%.3f",round(`Pr...z..`, digits = 3))))),
         z_value = paste0("z = ",paste0(sprintf("%.3f",round(`z.value`, digits = 3)))),
         stderr = paste0("(",as.character(format(round(`Boot.SE`, digits = 3),nsmall =3)),")"),
         Effectestimate = paste0(sprintf("%.3f",round(Estimate, digits = 3))," ",stderr),
         Effect = rownames(.)) %>% 
  select(Effect, Effectestimate, z_value, p_value) %>%
  pivot_longer(-c(Effect), names_to = "statistic", values_to = "value")

models_consent_boot <- full_join(model1,model2, by = c("Effect", "statistic")) %>% 
  mutate(n = c(49:51,1:6,16:48,7:15)) %>%
  #mutate(n = c(33,34,1:4,11:32,5:10)) 
  arrange(n) %>% 
  mutate(Effects = c(labels_boot)) %>% 
  select(Effects, model1 = value.x, model2 =value.y) %>% 
  rbind(c("N",as.character(nobs(mymodel_consent_total_con_inc)),as.character(nobs(mymodel_consent_factor_con_inc))),
        c("Log Likelihood", as.character(round(as.numeric(logLik(mymodel_consent_total_con_inc)),digits = 3)), as.character(round(as.numeric(logLik(mymodel_consent_factor_con_inc)),digits = 3))),
        c("AIC", as.character(round(mymodel_consent_total_con_inc$aic, digits = 3)), as.character(round(mymodel_consent_factor_con_inc$aic, digits = 3))))
```

```{r}
stargazer(models_consent_boot, type = "text", style = "default", align = F, summary=F, rownames = F, covariate.labels = c("","(1)","(2)"), table.placement = "H", no.space = T, single.row = F, header = F, dep.var.caption = NULL, title = "Consent")
```

## Table A21: Logistic Regression Models for the Outcome Variable Installation with Bootstrapped Standard Errors

```{r}
response <- "installation"
set.seed(123)

## conditional incentives binary

predictors_total_con_inc <- c("uncon_inc", "con_inc_total","pa_index","kh_index","age","sex","educ_medium","educ_high","married","hhsiz","german","rec_arm","study_tracking_num")

mymodel_installation_total_con_inc <- glm_installation_pa_kh_demo_coninc_dummy_consent


boot_results_installation_total_con_inc <- boot(data = gsurf_consentexp_analysis_full_consent, statistic = function(data, indices) boot_fn(gsurf_consentexp_analysis_full_consent, indices, response, predictors_total_con_inc), R = 1000)

# Extract bootstrapped standard errors
boot_se <- apply(boot_results_installation_total_con_inc$t, 2, sd)

# Display the original model coefficients with bootstrapped standard errors
original_coef <- coef(mymodel_installation_total_con_inc)
# Calculate z-values and p-values
z_values <- original_coef / boot_se
p_values <- round(2 * pnorm(-abs(z_values)), digits = 3)

# Create a results data frame
results_installation_total_con_inc <- data.frame(
  Estimate = original_coef,
  `Boot SE` = boot_se,
  `z value` = z_values,
  `Pr(>|z|)` = p_values
)

print(results_installation_total_con_inc)

## conditional incentives as factor
predictors_factor_con_inc <- c("uncon_inc", "as.factor(con_inc)","pa_index","kh_index","age","sex","educ_medium","educ_high","married","hhsiz","german","rec_arm","study_tracking_num")

mymodel_installation_factor_con_inc <-  glm_installation_pa_kh_demo_dummy_consent


boot_results_installation_factor_con_inc <- boot(data = gsurf_consentexp_analysis_full_consent, statistic = function(data, indices) boot_fn(gsurf_consentexp_analysis_full_consent, indices, response, predictors_factor_con_inc), R = 1000)

# Extract bootstrapped standard errors
boot_se <- apply(boot_results_installation_factor_con_inc$t, 2, sd)

# Display the original model coefficients with bootstrapped standard errors
original_coef <- coef(mymodel_installation_factor_con_inc)
# Calculate z-values and p-values
z_values <- original_coef / boot_se
p_values <- round(2 * pnorm(-abs(z_values)), digits = 3)

# Create a results data frame
results_installation_factor_con_inc <- data.frame(
  Estimate = original_coef,
  `Boot SE` = boot_se,
  `z value` = z_values,
  `Pr(>|z|)` = p_values
)

#print(results_installation_factor_con_inc)
```

```{r, include=FALSE}
labels_boot <- c("Unconditional 5€","","","Conditional any amount","","","Conditional 10€","","","Conditional 25€","","","Conditional 40€","","",
            "Privacy attitudes","","","Computer know-how","","","Age","","","Male","","","Education: medium (ref.: low)","","",
            "Education: high (ref.: low)","","","Married","","","Household size","", "","German nationality","","",
            "Recruitment via Meta (ref.: ALLBUS)","","","Study number (ref.: study two)","","", "Intercept","","")

model1 <- results_installation_total_con_inc %>% 
  mutate(p_value = if_else(`Pr...z..` <0.001, "p < 0.001",
                           paste0("p = ",paste0(sprintf("%.3f",round(`Pr...z..`, digits = 3))))),
         z_value = paste0("z = ",paste0(sprintf("%.3f",round(z.value, digits = 3)))),
         stderr = paste0("(",as.character(format(round(`Boot.SE`, digits = 3),nsmall =3)),")"),
         Effectestimate = paste0(sprintf("%.3f",round(Estimate, digits = 3))," ",stderr),
         Effect = rownames(.)) %>% 
  select(Effect, Effectestimate,z_value, p_value) %>%
  pivot_longer(-c(Effect), names_to = "statistic", values_to = "value")


model2 <- results_installation_factor_con_inc %>% 
  mutate(p_value = if_else(`Pr...z..` <0.001, "p < 0.001",
                           paste0("p = ",paste0(sprintf("%.3f",round(`Pr...z..`, digits = 3))))),
         z_value = paste0("z = ",paste0(sprintf("%.3f",round(z.value, digits = 3)))),
         stderr = paste0("(",as.character(format(round(`Boot.SE`, digits = 3),nsmall =3)),")"),
         Effectestimate = paste0(sprintf("%.3f",round(Estimate, digits = 3))," ",stderr),
         Effect = rownames(.)) %>% 
  select(Effect, Effectestimate,z_value, p_value) %>%
  pivot_longer(-c(Effect), names_to = "statistic", values_to = "value")

models_installation_boot <- full_join(model1,model2, by = c("Effect", "statistic")) %>% 
  mutate(n = c(49:51,1:6,16:48,7:15)) %>% 
  arrange(n) %>% 
  mutate(Effects = c(labels_boot)) %>% 
  select(Effects, model1 = value.x, model2 =value.y) %>% 
  rbind(c("N",as.character(nobs(mymodel_installation_total_con_inc)),as.character(nobs(mymodel_installation_factor_con_inc))),
        c("Log Likelihood", as.character(round(as.numeric(logLik(mymodel_installation_total_con_inc)),digits = 3)), as.character(round(as.numeric(logLik(mymodel_installation_factor_con_inc)),digits = 3))),
        c("AIC", as.character(round(mymodel_installation_total_con_inc$aic, digits = 3)), as.character(round(mymodel_installation_factor_con_inc$aic, digits = 3))))
```

```{r}
stargazer(models_installation_boot, type = "text", style = "default", align = F, summary=F, rownames = F, covariate.labels = c("","(1)","(2)"), table.placement = "H", no.space = T, single.row = F, header = F, dep.var.caption = NULL, title = "Installation")
```

## Table A22: Negative Binomial Regression with Bootstrapped Standard Errors

```{r}
# Define a function for bootstrapping
boot_fn_glmmTMB <- function(data, i) {
  # Bootstrap the data
  d <- data[i, ]
  # Fit the glmmTMB model
  model <- glmmTMB(active_days_threshold  ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, 
                   data = d, family = nbinom2)
  # Extract the conditional model coefficients
  cf <- unlist(fixef(model))
  return(cf)
}

```

```{r}
# conditional incentive total (binary)
mymodel_active_days_total_coninc <- lm_active_full_transfer_coninc_dummy 

boot_results_active_days_total_coninc <- boot(data = gsurf_consentexp_analysis_full_active, statistic = boot_fn_glmmTMB, R = 1000)

# Extract bootstrapped standard errors
boot_se <- apply(boot_results_active_days_total_coninc$t, 2, sd)

# Display the original model coefficients with bootstrapped standard errors
original_coef <- as.data.frame(summary(mymodel_active_days_total_coninc)$coefficients$cond)$Estimate
# Calculate z-values and p-values
z_values <- original_coef / boot_se[1:14]
p_values <- round(2 * pnorm(-abs(z_values)), digits = 3)

# Create a results data frame
results_active_days_total_coninc <- data.frame(
  Estimate = original_coef,
  `Boot SE` = boot_se[1:14],
  `z value` = z_values,
  `Pr(>|z|)` = p_values
)
rownames(results_active_days_total_coninc) <- names(results_active_days_total_coninc$t0)[1:16]

print(results_active_days_total_coninc)

# conditional incentive as factor

mymodel_active_days_factor_coninc <- lm_active_full_transfer_demo 
boot_results_active_days_factor_coninc <- boot(data = gsurf_consentexp_analysis_full_active, statistic = boot_fn_glmmTMB, R = 1000)

# Extract bootstrapped standard errors
boot_se <- apply(boot_results_active_days_factor_coninc$t, 2, sd)

# Display the original model coefficients with bootstrapped standard errors
original_coef <- as.data.frame(summary(mymodel_active_days_factor_coninc)$coefficients$cond)$Estimate
# Calculate z-values and p-values
z_values <- original_coef / boot_se[1:16]
p_values <- round(2 * pnorm(-abs(z_values)), digits = 3)

# Create a results data frame
results_active_days_factor_coninc <- data.frame(
  Estimate = original_coef,
  `Boot SE` = boot_se[1:16],
  `z value` = z_values,
  `Pr(>|z|)` = p_values
)
rownames(results_active_days_factor_coninc) <- names(results_active_days_factor_coninc$t0)[1:16]

#print(results_active_days_factor_coninc)
```

```{r, include=FALSE}
labels_boot <- c("Unconditional 5€","","","Conditional any amount","","","Conditional 10€","","","Conditional 25€","","","Conditional 40€","","",
            "Privacy attitudes","","","Computer know-how","","","Age","","","Male","","","Education: medium (ref.: low)","","",
            "Education: high (ref.: low)","","","Married","","","Household size","", "","German nationality","","",
            "Recruitment via Meta (ref.: ALLBUS)","","","Study number (ref.: study two)","","", "Intercept","","")

model1 <- results_active_days_total_coninc %>% 
  mutate(p_value = if_else(`Pr...z..` <0.001, "p < 0.001",
                           paste0("p = ",paste0(sprintf("%.3f",round(`Pr...z..`, digits = 3))))),
         z_value = paste0("z = ",paste0(sprintf("%.3f",round(z.value, digits = 3)))),
         stderr = paste0("(",as.character(format(round(`Boot.SE`, digits = 3),nsmall =3)),")"),
         Effectestimate = paste0(sprintf("%.3f",round(Estimate, digits = 3))," ",stderr),
         Effect = c("Intercept", "uncon_inc", "con_inc_total","pa_index","kh_index","age","sex","educ_medium","educ_high","married","hhsiz","german","rec_arm","study_tracking_num")) %>% 
  select(Effect, Effectestimate,z_value, p_value) %>%
  pivot_longer(c(Effectestimate, z_value, p_value), names_to = "statistic", values_to = "value")


model2 <- results_active_days_factor_coninc %>% 
  mutate(p_value = if_else(`Pr...z..` <0.001, "p < 0.001",
                           paste0("p = ",paste0(sprintf("%.3f",round(`Pr...z..`, digits = 3))))),
         z_value = paste0("z = ",paste0(sprintf("%.3f",round(z.value, digits = 3)))),
         stderr = paste0("(",as.character(format(round(`Boot.SE`, digits = 3),nsmall =3)),")"),
         Effectestimate = paste0(sprintf("%.3f",round(Estimate, digits = 3))," ",stderr),
         Effect = c("Intercept", "uncon_inc", "con_inc10","con_inc25","con_inc40","pa_index","kh_index","age","sex","educ_medium","educ_high","married","hhsiz","german","rec_arm","study_tracking_num")) %>% 
  select(Effect, Effectestimate,z_value, p_value) %>%
  pivot_longer(c(Effectestimate,z_value, p_value), names_to = "statistic", values_to = "value")


models_tracked_days_boot <- full_join(model1,model2, by = c("Effect", "statistic")) %>% 
  mutate(n = c(49:51,1:6,16:48,7:15)) %>% 
  arrange(n) %>% 
  mutate(Effects = c(labels_boot)) %>% 
  select(Effects, model1 = value.x, model2 =value.y) %>%
  rbind(c("N",as.character(nobs(mymodel_active_days_total_coninc)),as.character(nobs(mymodel_active_days_factor_coninc))),
        c("Log Likelihood", as.character(round(summary(mymodel_active_days_total_coninc)$AIC[3],digits = 3)), as.character(round(summary(mymodel_active_days_factor_coninc)$AIC[3],digits = 3))),
        c("AIC", as.character(round(summary(mymodel_active_days_total_coninc)$AIC[1],digits = 3)), as.character(round(summary(mymodel_active_days_factor_coninc)$AIC[1],digits = 3))))

```

```{r}
stargazer(models_tracked_days_boot, type = "text", style = "default", align = F, summary=F, rownames = F, covariate.labels = c("","(1)","(2)"), table.placement = "H", no.space = T, single.row = F, header = F, dep.var.caption = NULL, title = "Tracked Days")
```

## Table A23: Logistic Regression Models Using Backward Difference Coding

```{r}
glm_consent_pa_kh_demo_backward <- glm(consent ~ uncon_inc+con_inc_backward+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full, family = "binomial")

glm_installation_pa_kh_demo_backward <- glm(installation ~ uncon_inc+con_inc_backward+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full_consent, family = "binomial")
```

```{r}
labels_backward <- c("Unconditional 5€","Conditional 10€ (ref.: 0€)","Conditional 25€ (ref.: 10€)","Conditional 40€ (ref.: 25€",
            "Privacy attitudes","Computer know-how","Age","Male","Education: medium (ref.: low)",
            "Education: high (ref.: low)","Married","Household size", "German nationality",
            "Recruitment via Meta (ref.: ALLBUS)","Study number (ref.: study two)", "Intercept")

stargazer(glm_consent_pa_kh_demo_backward, glm_installation_pa_kh_demo_backward, type = "text", style = "default", align = F, summary=T, rownames = T, covariate.labels = labels_backward, table.placement = "H", no.space = T, single.row = T, header = F, star.cutoffs = NA, dep.var.caption = NULL, report = "vcstp", dep.var.labels = c("Consent","Installation"))
```


## Table A24: Negative Binomial Regression Using Backward Difference Coding

```{r}
lm_active_full_transfer_demo_backward <- glmmTMB(active_days_threshold  ~ uncon_inc+con_inc_backward+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, family=nbinom2, data = gsurf_consentexp_analysis_full_active)
```

```{r, include=FALSE}
labels_backward <- c("Unconditional 5€","", "","Conditional 10€ (ref.: 0€)","","","Conditional 25€ (ref.: 10€)","","","Conditional 40€ (ref.: 25€","","",
                     "Privacy attitudes","","","Computer know-how","","","Age","","","Male","","","Education: medium (ref.: low)","","",
            "Education: high (ref.: low)","","","Married","","","Household size","","", "German nationality","","",
            "Recruitment via Meta (ref.: ALLBUS)","","","Study number (ref.: study two)","","","Intercept","","")

model1a <- data.frame(Estimate = unlist(lapply(coef(summary(lm_active_full_transfer_demo_backward)),function(x) x[,"Estimate"])),
                      stderr = paste0("(",as.character(round(unlist(lapply(coef(summary(lm_active_full_transfer_demo_backward)),function(x) x[,"Std. Error"])),3)),")"),
                      zvalue = unlist(lapply(coef(summary((lm_active_full_transfer_demo_backward))),function(x) x[,"z value"])),
                      pvalue = unlist(lapply(coef(summary(lm_active_full_transfer_demo_backward)),function(x) x[,"Pr(>|z|)"]))) %>%
  mutate(Effect = rownames(.),
         p_value = if_else(pvalue <0.001, "p < 0.001",
                           paste0("p = ",paste0(sprintf("%.3f",round(pvalue, digits = 3))))),
         z_value = paste0("z = ",paste0(sprintf("%.3f",round(zvalue, digits = 3)))),
         Effectestimate = paste0(sprintf("%.3f",round(Estimate, digits = 3))," ",stderr)) %>% 
  select(Effect, Effectestimate,z_value, p_value) %>% 
  pivot_longer(-c(Effect), names_to = "statistic", values_to = "value")


models_active_days_backward <- model1a %>% 
  mutate(n = c(46:48,1:45)) %>% 
  arrange(n) %>% 
  mutate(Effects = c(labels_backward)) %>% 
  select(Effects, model1 = value) %>% 
  rbind(c("N",as.character(nobs(lm_active_full_transfer_demo_backward))),
        c("Log Likelihood", as.character(round(summary(lm_active_full_transfer_demo_backward)$AIC[3],digits = 3))),
        c("AIC", as.character(round(summary(lm_active_full_transfer_demo_backward)$AIC[1],digits = 3))))


```

```{r}
stargazer(models_active_days_backward, type = "text", style = "default", align = F, summary=F, rownames = F, covariate.labels = c("","(1)","(2)"), table.placement = "H", no.space = T, single.row = F, header = F, dep.var.caption = NULL, title = "Tracked Days")
```

## Table A25: Supplementary Analysis: Logistic Regression for Installation (Compound Outcome: Consent + Installation; Full Sample)

```{r}
gsurf_consentexp_analysis_full <- gsurf_consentexp_analysis_full %>% 
  mutate(installation_total = ifelse(is.na(installation), 0, installation))

glm_installation_pa_kh_demo_coninc_dummy <- glm(installation_total ~ uncon_inc+con_inc_total+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full, family = "binomial")

glm_installation_pa_kh_demo_dummy <- glm(installation_total ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full, family = "binomial")

labels <- c("Unconditional 5€","Conditional any amount","Conditional 10€","Conditional 25€","Conditional 40€",
            "Privacy attitudes","Computer know-how","Age","Male","Education: medium (ref.: low)",
            "Education: high (ref.: low)","Married","Household size", "German nationality",
            "Recruitment via Meta (ref.: ALLBUS)","Study number (ref.: study two)", "Intercept")

stargazer(glm_installation_pa_kh_demo_coninc_dummy, glm_installation_pa_kh_demo_dummy, type = "text", style = "default", align = F, summary=T, rownames = T, covariate.labels = labels, table.placement = "H", no.space = T, single.row = T, header = F, star.cutoffs = NA, dep.var.caption = NULL, report = "vcstp", dep.var.labels = "Installation")
```


## Table A26: Supplementary Analysis: Logistic Regression for Sustained Participation (Compound Outcome: Consent + Installation + ≥ 30 Active Days; Full Sample)

```{r}
glm_30day_pa_kh_demo_coninc_dummy <- glm(threshold30days ~ uncon_inc+con_inc_total+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full, family = "binomial")

glm_30day_pa_kh_demo_dummy <- glm(threshold30days ~ uncon_inc+as.factor(con_inc)+pa_index+kh_index+age+sex+educ_medium+educ_high+married+hhsiz+german+rec_arm+study_tracking_num, data = gsurf_consentexp_analysis_full, family = "binomial")

labels <- c("Unconditional 5€","Conditional any amount","Conditional 10€","Conditional 25€","Conditional 40€",
            "Privacy attitudes","Computer know-how","Age","Male","Education: medium (ref.: low)",
            "Education: high (ref.: low)","Married","Household size", "German nationality",
            "Recruitment via Meta (ref.: ALLBUS)","Study number (ref.: study two)", "Intercept")

stargazer(glm_30day_pa_kh_demo_coninc_dummy, glm_30day_pa_kh_demo_dummy, type = "text", style = "default", align = F, summary=T, rownames = T, covariate.labels = labels, table.placement = "H", no.space = T, single.row = T, header = F, star.cutoffs = NA, dep.var.caption = NULL, report = "vcstp", dep.var.labels = "Threshold 30 tracked days reached")
```


# Supplementrary analysis of Cost Efficiency

```{r, warning=FALSE}
#### installation (full sample)

## Unconditional incentive

# Get estimated marginal means (EMMs)
emm_installation_uncon_inc_full_sample <- emmeans(glm_installation_pa_kh_demo_coninc_dummy, ~ uncon_inc, type = "response", weights = "proportional")
summary(emm_installation_uncon_inc_full_sample, level = 0.834) #for CIs at 83.4% intervalls

## Conditional incentive

# Get estimated marginal means (EMMs)
emm_installation_con_inc_full_sample <- emmeans(glm_installation_pa_kh_demo_dummy, ~ con_inc, type = "response", weights = "proportional")
summary(emm_installation_con_inc_full_sample, level = 0.834) #for CIs at 83.4% intervalls

#### tracked days (full sample)

## Unconditional incentive

# Get estimated marginal means (EMMs)
emm_30day_uncon_inc_full_sample <- emmeans(glm_30day_pa_kh_demo_coninc_dummy, ~ uncon_inc, type = "response", weights = "proportional")
summary(emm_30day_uncon_inc_full_sample, level = 0.834) #for CIs at 83.4% intervalls

## Conditional incentive

# Get estimated marginal means (EMMs)
emm_30day_con_inc_full_sample <- emmeans(glm_30day_pa_kh_demo_dummy, ~ con_inc, type = "response", weights = "proportional")
summary(emm_30day_con_inc_full_sample, level = 0.834) #for CIs at 83.4% intervalls
summary(emm_30day_con_inc_full_sample, level = 0.834)$prob #for exact values

```
## Table A27: Estimated Direct Costs of Incentivizing 1,000 Participants by Incentive Condition
```{r}
cost5_1000 <- 5*1000
cost10_1000 <- summary(emm_30day_con_inc_full_sample)$prob[2]*10*1000
cost25_1000 <- summary(emm_30day_con_inc_full_sample)$prob[3]*25*1000
cost40_1000 <- summary(emm_30day_con_inc_full_sample)$prob[4]*40*1000


cost_assessment <- data_frame(amount = c("5€","10€","25€","40€"),
                              condition = c("Unconditional","Conditional","Conditional","Conditional"),
                              icentives_1000_participants_in_euro = round(c(cost5_1000,cost10_1000,cost25_1000,cost40_1000),0))
cost_assessment
```

## Results Reported in Text

financial cost of achieving 100 additional participants at different participation stages
```{r}
#consent: drawing on probabilities estimated on model presented in table A9
cost5_100_consent <- (cost5_1000 / (summary(emm_consent_uncon_inc)$prob[2]*1000 - summary(emm_consent_uncon_inc)$prob[1]*1000))*100
cost10_100_consent <- (cost10_1000 / (summary(emm_consent_con_inc)$prob[2]*1000 - summary(emm_consent_con_inc)$prob[1]*1000))*100
cost25_100_consent <- (cost25_1000 / (summary(emm_consent_con_inc)$prob[3]*1000 - summary(emm_consent_con_inc)$prob[1]*1000))*100
cost40_100_consent <- (cost40_1000 / (summary(emm_consent_con_inc)$prob[4]*1000 - summary(emm_consent_con_inc)$prob[1]*1000))*100

#installation in the complete cohort: drawing on probabilities estimated on model presented in table A22
cost5_100_installation <- (cost5_1000 / (summary(emm_installation_uncon_inc_full_sample)$prob[2]*1000 - summary(emm_installation_uncon_inc_full_sample)$prob[1]*1000))*100
cost10_100_installation <- (cost10_1000 / (summary(emm_installation_con_inc_full_sample)$prob[2]*1000 - summary(emm_installation_con_inc_full_sample)$prob[1]*1000))*100
cost25_100_installation <- (cost25_1000 / (summary(emm_installation_con_inc_full_sample)$prob[3]*1000 - summary(emm_installation_con_inc_full_sample)$prob[1]*1000))*100
cost40_100_installation <- (cost40_1000 / (summary(emm_installation_con_inc_full_sample)$prob[4]*1000 - summary(emm_installation_con_inc_full_sample)$prob[1]*1000))*100

#installation those who gave consent as presented in table A10
cost5_100_installation_consent <- (cost5_1000 / (summary(emm_installation_uncon_inc)$prob[2]*1000 - summary(emm_installation_uncon_inc)$prob[1]*1000))*100
cost10_100_installation_consent <- (cost10_1000 / (summary(emm_installation_con_inc)$prob[2]*1000 - summary(emm_installation_con_inc)$prob[1]*1000))*100
cost25_100_installation_consent <- (cost25_1000 / (summary(emm_installation_con_inc)$prob[3]*1000 - summary(emm_installation_con_inc)$prob[1]*1000))*100
cost40_100_installation_consent <- (cost40_1000 / (summary(emm_installation_con_inc)$prob[4]*1000 - summary(emm_installation_con_inc)$prob[1]*1000))*100


#30 tracked days (sustained engagement) in the complete cohort: drawing on probabilities estimated on model presented in table A23
cost5_100_30days <- (cost5_1000 / (summary(emm_30day_uncon_inc_full_sample)$prob[2]*1000 - summary(emm_30day_uncon_inc_full_sample)$prob[1]*1000))*100
cost10_100_30days <- (cost10_1000 / (summary(emm_30day_con_inc_full_sample)$prob[2]*1000 - summary(emm_30day_con_inc_full_sample)$prob[1]*1000))*100
cost25_100_30days <- (cost25_1000 / (summary(emm_30day_con_inc_full_sample)$prob[3]*1000 - summary(emm_30day_con_inc_full_sample)$prob[1]*1000))*100
cost40_100_30days <- (cost40_1000 / (summary(emm_30day_con_inc_full_sample)$prob[4]*1000 - summary(emm_30day_con_inc_full_sample)$prob[1]*1000))*100

additional_outcomes100 <- data.frame(amount = c("5€","10€","25€","40€"),
                                     condition = c("Unconditional","Conditional","Conditional","Conditional"),
                                     consent_100additional_in_euro = c(cost5_100_consent,cost10_100_consent,cost25_100_consent,cost40_100_consent),
                                     installation_100additional_in_euro = c(cost5_100_installation,cost10_100_installation,cost25_100_installation,cost40_100_installation),
                                     installation_consent_100additional_in_euro = c(cost5_100_installation_consent,cost10_100_installation_consent,
                                                                                    cost25_100_installation_consent,cost40_100_installation_consent),
                                     days30_100additional_in_euro = c(cost5_100_30days,cost10_100_30days,cost25_100_30days,cost40_100_30days))
additional_outcomes100
```

the number of additional outcomes achievable per €10,000 invested
```{r}
#consent: drawing on probabilities estimated on model presented in table A9
gain5_consent <- round(10000 / (cost5_100_consent/100),0)
gain10_consent <- round(10000 / (cost10_100_consent/100),0)
gain25_consent <- round(10000 / (cost25_100_consent/100),0)
gain40_consent <- round(10000 / (cost40_100_consent/100),0)

#installation in the complete cohort: drawing on probabilities estimated on model presented in table A22
gain5_installation <- round(10000 / (cost5_100_installation/100),0)
gain10_installation <- round(10000 / (cost10_100_installation/100),0)
gain25_installation <- round(10000 / (cost25_100_installation/100),0)
gain40_installation <- round(10000 / (cost40_100_installation/100),0)

#installation those who gave consent as presented in Table A10
gain5_installation_consent <- round(10000 / (cost5_100_installation_consent/100),0)
gain10_installation_consent <- round(10000 / (cost10_100_installation_consent/100),0)
gain25_installation_consent <- round(10000 / (cost25_100_installation_consent/100),0)
gain40_installation_consent <- round(10000 / (cost40_100_installation_consent/100),0)


#30 tracked days (sustained engagement) in the complete cohort: drawing on probabilities estimated on model presented in table A23
gain5_30days <- round(10000 / (cost5_100_30days/100),0)
gain10_30days <- round(10000 / (cost10_100_30days/100),0)
gain25_30days <- round(10000 / (cost25_100_30days/100),0)
gain40_30days <- round(10000 / (cost40_100_30days/100),0)

gain10000euros <- data.frame(amount = c("5€","10€","25€","40€"),
                                     condition = c("Unconditional","Conditional","Conditional","Conditional"),
                                     consent_100additional_in_euro = c(gain5_consent,gain10_consent,gain25_consent,gain40_consent),
                                     installation_100additional_in_euro = c(gain5_installation,gain10_installation,
                                                                            gain25_installation,gain40_installation),
                             installation_consent_100additional_in_euro = c(gain5_installation_consent,gain10_installation_consent,
                                                                            gain25_installation_consent,gain40_installation_consent),
                                     days30_100additional_in_euro = c(gain5_30days,gain10_30days,gain25_30days,gain40_30days))
gain10000euros
```


# Figures

```{r, echo=FALSE}
df_emm_consent_F1 <-  data.frame(outcome = "consent",
                               incentive = c("unconditional","unconditional","conditional","conditional"),
                               promised = c(0,1,0,1), #"0","5","0","10","25","40"
                               ame_predict =c(summary(emm_consent_uncon_inc, level = 0.834)$prob,summary(emm_consent_con_inc_dummy, level = 0.834)$prob),
                               conf.low = c(summary(emm_consent_uncon_inc, level = 0.834)$asymp.LCL,summary(emm_consent_con_inc_dummy, level = 0.834)$asymp.LCL),
                               conf.high = c(summary(emm_consent_uncon_inc, level = 0.834)$asymp.UCL,summary(emm_consent_con_inc_dummy, level = 0.834)$asymp.UCL))

df_emm_installation_F1 <- data.frame(outcome = "installation",
                               incentive = c("unconditional","unconditional","conditional","conditional"),
                               promised = c(0,1,0,1),
                               ame_predict =c(summary(emm_installation_uncon_inc, level = 0.834)$prob,summary(emm_installation_con_inc_dummy, level = 0.834)$prob),
                               conf.low = c(summary(emm_installation_uncon_inc, level = 0.834)$asymp.LCL,summary(emm_installation_con_inc_dummy, level = 0.834)$asymp.LCL),
                               conf.high = c(summary(emm_installation_uncon_inc, level = 0.834)$asymp.UCL,summary(emm_installation_con_inc_dummy, level = 0.834)$asymp.UCL))

df_emm_active_F1 <- data.frame(outcome = "active days",
                               incentive = c("unconditional","unconditional","conditional","conditional"),
                               promised = c(0,1,0,1),
                               ame_predict =c(summary(emm_activdays_uncon_inc, level = 0.834)$response,summary(emm_activdays_con_inc_dummy, level = 0.834)$response),
                               conf.low = c(summary(emm_activdays_uncon_inc, level = 0.834)$asymp.LCL,summary(emm_activdays_con_inc_dummy, level = 0.834)$asymp.LCL),
                               conf.high = c(summary(emm_activdays_uncon_inc, level = 0.834)$asymp.UCL,summary(emm_activdays_con_inc_dummy, level = 0.834)$asymp.UCL))

df_emmF1 <- rbind(df_emm_consent_F1,df_emm_installation_F1,df_emm_active_F1) %>% 
  mutate(type =if_else(promised == 0,"no incentive","incentive")) %>% 
  mutate(type = factor(type, levels=unique(type))) %>% 
  mutate(Zero = if_else(promised == 0,ame_predict,NA)) %>% 
  fill(Zero) %>% 
  mutate(Zero = if_else(promised == 0,NA,Zero)) %>% 
  mutate(conf.low = round(conf.low, digits = 3),
         conf.high = round(conf.high, digits = 3),
         ame_predict = round(ame_predict, digits = 3))
```

```{r}
scaleFUN <- function(x){if_else(x < 1,sprintf("%.2f", x),sprintf("%.0f", x))} 

df_emmF1c <- df_emmF1 %>% 
  filter(outcome == "consent") %>% 
  mutate(incentive = case_match(incentive,
                                "unconditional" ~ "Unconditional",
                                "conditional" ~ "Conditional"))
p4 <-
ggplot(df_emmF1c, aes(x = type, y = ame_predict, color = as.factor(type))) +
  geom_segment(x = 1,xend = 2,y=df_emmF1c$Zero,yend=df_emmF1c$ame_predict, data = df_emmF1c, size = 2, show.legend=F) + #
  geom_point(size = 5) +
  geom_errorbar(data=df_emmF1c,aes(ymin=conf.low,ymax=conf.high),
                width = 0.4, size = 1,
                show.legend = FALSE)+
  ylab("Predicted Probability of Consent")+
  xlab("")+
  scale_y_continuous(labels=scaleFUN)+
  scale_x_discrete(labels = c("No\nIncentive","Incentive\n"))+
  scale_color_manual(name = "Incentive", labels = c("no","yes"), values = c("#000000", "#000000"))+
  ggh4x::facet_grid2(factor(incentive, levels = c("Unconditional","Conditional")) ~ factor(outcome, levels = c("consent", "installation", "active days"),
                                                                                           labels = c("Consent", "Installation", "tracked days")), scales = "free", independent = "y", switch = "y")+
  ggh4x::scale_y_facet(COL == 1, limits = c(0.35,0.60), breaks = c(.4,.5,.6), labels = scales::percent) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none")

df_emmF1i <- df_emmF1 %>% 
  filter(outcome == "installation")%>% 
  mutate(incentive = case_match(incentive,
                                "unconditional" ~ "Unconditional",
                                "conditional" ~ "Conditional"))
p5 <-
  ggplot(df_emmF1i, aes(x = type, y = ame_predict, color = as.factor(type))) +
  geom_segment(x = 1,xend = 2,y=df_emmF1i$Zero,yend=df_emmF1i$ame_predict, data = df_emmF1i, size = 2, show.legend=F) + #
  geom_point(size = 5) +
  geom_errorbar(data=df_emmF1i,aes(ymin=conf.low,ymax=conf.high),
                width = 0.4, size = 1,
                show.legend = FALSE)+
  ylab("Predicted Probability of Installation")+
  xlab("")+
  scale_y_continuous(labels=scaleFUN)+
  scale_x_discrete(labels = c("No\nIncentive","Incentive\n"))+
  scale_color_manual(name = "Incentive", labels = c("no","yes"), values = c("#000000", "#000000"))+
  ggh4x::facet_grid2(factor(incentive, levels = c("Unconditional","Conditional")) ~ factor(outcome, levels = c("consent", "installation", "active days"), labels = c("Consent", "Installation", "tracked days")), scales = "free", independent = "y", switch = "y")+
  ggh4x::scale_y_facet(COL == 1, limits = c(0.35,0.60), breaks = c(.4,.5,.6), labels = scales::percent) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none")


df_emmF1a <- df_emmF1 %>% 
  filter(outcome == "active days") %>% 
  mutate(incentive = case_match(incentive,
                                "unconditional" ~ "Unconditional",
                                "conditional" ~ "Conditional"))
p6 <-
  ggplot(df_emmF1a, aes(x = type, y = ame_predict, color = as.factor(type))) +
  geom_segment(x = 1,xend = 2,y=df_emmF1a$Zero,yend=df_emmF1a$ame_predict, data = df_emmF1a, size = 2, show.legend=F) + #
  geom_point(size = 5) +
  geom_errorbar(data=df_emmF1a,aes(ymin=conf.low,ymax=conf.high),
                width = 0.4, size = 1,
                show.legend = FALSE)+
  ylab("Predicted Number of Tracked Days")+
  xlab("")+
  scale_y_continuous(labels=scaleFUN)+
  scale_x_discrete(labels = c("No\nIncentive","Incentive\n"))+
  scale_color_manual(name = "Incentive", labels = c("no","yes"), values = c("#000000", "#000000"))+
  ggh4x::facet_grid2(factor(incentive, levels = c("Unconditional","Conditional")) ~ factor(outcome, levels = c("consent", "installation", "active days"),
                                                                                           labels = c("consent", "installation", "Number of Tracked Days")),
                     scales = "free", independent = "y", switch = "y")+
  ggh4x::scale_y_facet(COL == 1, limits = c(24,40)) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none")

grid.arrange(p4,p5,p6, nrow = 1)
```

```{r, echo=FALSE}
df_emm_consent_F2 <-  data.frame(outcome = "consent",
                               incentive = c("conditional","conditional","conditional","conditional"),
                               amount = c(0,10,25,40), #"0","5","0","10","25","40"
                               ame_predict =c(summary(emm_consent_con_inc, level = 0.834)$prob),
                               conf.low = c(summary(emm_consent_con_inc, level = 0.834)$asymp.LCL),
                               conf.high = c(summary(emm_consent_con_inc, level = 0.834)$asymp.UCL))

df_emm_installation_F2 <- data.frame(outcome = "installation",
                               incentive = c("conditional","conditional","conditional","conditional"),
                               amount = c(0,10,25,40),
                               ame_predict =c(summary(emm_installation_con_inc, level = 0.834)$prob),
                               conf.low = c(summary(emm_installation_con_inc, level = 0.834)$asymp.LCL),
                               conf.high = c(summary(emm_installation_con_inc, level = 0.834)$asymp.UCL))

df_emm_active_F2 <- data.frame(outcome = "active days",
                               incentive = c("conditional","conditional","conditional","conditional"),
                               amount = c(0,10,25,40),
                               ame_predict =c(summary(emm_activdays_con_inc, level = 0.834)$response),
                               conf.low = c(summary(emm_activdays_con_inc, level = 0.834)$asymp.LCL),
                               conf.high = c(summary(emm_activdays_con_inc, level = 0.834)$asymp.UCL))

df_emm_F2 <- rbind(df_emm_consent_F2,df_emm_installation_F2,df_emm_active_F2) %>% 
  mutate(type =if_else(amount == 0,"no incentive",paste0(as.character(amount),"€ incentive"))) %>% 
  mutate(type = factor(type, levels=unique(type))) %>% 
  mutate(Zero = if_else(amount == 0,ame_predict,NA)) %>% 
  fill(Zero) %>% 
  mutate(Zero = if_else(amount == 0,NA,Zero)) %>% 
  mutate(conf.low = round(conf.low, digits = 3),
         conf.high = round(conf.high, digits = 3),
         ame_predict = round(ame_predict, digits = 3))
```

```{r, echo=FALSE}
df_emm_F3 <- df_emm_F2 %>% 
  filter(outcome == "consent") %>% 
  mutate(outcome = "Consent")

p1 <-
ggplot(df_emm_F3, aes(x = as.factor(amount), y = ame_predict, color = as.factor(amount))) +
  geom_segment(x = 1,xend = as.factor(df_emm_F3$amount),y=df_emm_F3$Zero,yend=df_emm_F3$ame_predict, data = df_emm_F3, size = 2, show.legend=F) +
  geom_point(size = 5) +
  geom_errorbar(data=df_emm_F3,aes(ymin=conf.low,ymax=conf.high),
                width = 0.4, size = 1,
                show.legend = FALSE)+
  ylab("Predicted Probability of Consent")+
  xlab("Incentive Height")+
  scale_y_continuous(labels=scaleFUN)+
  scale_color_manual(name = "Incentive amount", labels = c("0€","10€","25€","40€"), values = c("#000000", "#000000",
                                                                                               "#000000", "#000000"))+
  scale_x_discrete(labels = c("0€","10€","25€","40€")) +
  ggh4x::facet_grid2(factor(incentive,
                            levels = c("conditional")) ~ outcome,
                     scales = "free", independent = "y", switch = "y", remove_labels = "x")+
  ggh4x::scale_y_facet(COL == 1, limits = c(0.35,0.70), labels = scales::percent) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none")

df_emm_F3 <- df_emm_F2 %>% 
  filter(outcome == "installation") %>% 
  mutate(outcome = "Installation")

p2 <-
  ggplot(df_emm_F3, aes(x = as.factor(amount), y = ame_predict, color = as.factor(amount))) +
  geom_segment(x = 1,xend = as.factor(df_emm_F3$amount),y=df_emm_F3$Zero,yend=df_emm_F3$ame_predict, data = df_emm_F3, size = 2, show.legend=F) +
  geom_point(size = 5) +
  geom_errorbar(data=df_emm_F3,aes(ymin=conf.low,ymax=conf.high),
                width = 0.4, size = 1,
                show.legend = FALSE)+
  ylab("Predicted Probability of Installation")+
  xlab("Incentive Height")+
  scale_y_continuous(labels=scaleFUN)+
  scale_color_manual(name = "Incentive amount", labels = c("0€","10€","25€","40€"), values = c("#000000", "#000000",
                                                                                               "#000000", "#000000"))+
  scale_x_discrete(labels = c("0€","10€","25€","40€")) +
  ggh4x::facet_grid2(factor(incentive,
                            levels = c("conditional")) ~ outcome, 
                     scales = "free", independent = "y", switch = "y", remove_labels = "x")+
  ggh4x::scale_y_facet(COL == 1, limits = c(0.35,0.70), labels = scales::percent) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none")

df_emm_F3 <- df_emm_F2 %>% 
  filter(outcome == "active days")

p3 <-
  ggplot(df_emm_F3, aes(x = as.factor(amount), y = ame_predict, color = as.factor(amount))) +
  geom_segment(x = 1,xend = as.factor(df_emm_F3$amount),y=df_emm_F3$Zero,yend=df_emm_F3$ame_predict, data = df_emm_F3, size = 2, show.legend=F) +
  geom_point(size = 5) +
  geom_errorbar(data=df_emm_F3,aes(ymin=conf.low,ymax=conf.high),
                width = 0.4, size = 1,
                show.legend = FALSE)+
  ylab("Predicted Number of Tracked Days")+
  xlab("Incentive Height")+
  scale_y_continuous(labels=scaleFUN)+
  scale_color_manual(name = "Incentive amount", labels = c("0€","10€","25€","40€"), values = c("#000000", "#000000",
                                                                                               "#000000", "#000000"))+
  scale_x_discrete(labels = c("0€","10€","25€","40€")) +
  ggh4x::facet_grid2(factor(incentive,
                            levels = c("conditional")) ~ factor(outcome, levels = "active days", labels = "Number of Tracked Days"),  
                     scales = "free", independent = "y", switch = "y", remove_labels = "x")+

  ggh4x::scale_y_facet(COL == 1, limits = c(24,40)) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none")

grid.arrange(p1,p2,p3, nrow = 1)

```


```{r}
df_emm_consent_F4 <- data.frame(outcome = "Consent",
                               recruitment_arm = c("ALLBUS","Meta"),
                               ame_predict =c(summary(emm_consent_rec_arm, level = 0.834)$prob),
                               conf.low = c(summary(emm_consent_rec_arm, level = 0.834)$asymp.LCL),
                               conf.high = c(summary(emm_consent_rec_arm, level = 0.834)$asymp.UCL))

df_emm_installation_F4 <- data.frame(outcome = "Installation",
                               recruitment_arm = c("ALLBUS","Meta"),
                               ame_predict =c(summary(emm_installation_rec_arm, level = 0.834)$prob),
                               conf.low = c(summary(emm_installation_rec_arm, level = 0.834)$asymp.LCL),
                               conf.high = c(summary(emm_installation_rec_arm, level = 0.834)$asymp.UCL))

df_emm_active_F4 <- data.frame(outcome = "Number of Tracked Days",
                               recruitment_arm = c("ALLBUS","Meta"),
                               ame_predict =c(summary(emm_activdays_rec_arm, level = 0.834)$response),
                               conf.low = c(summary(emm_activdays_rec_arm, level = 0.834)$asymp.LCL),
                               conf.high = c(summary(emm_activdays_rec_arm, level = 0.834)$asymp.UCL))
```

```{r}
px <-
ggplot(df_emm_consent_F4, aes(x = rev(as.factor(recruitment_arm)), y = rev(ame_predict), color = rev(as.factor(recruitment_arm)))) +
  geom_point(size = 5) +
  geom_errorbar(data=df_emm_consent_F4,aes(ymin=rev(conf.low),ymax=rev(conf.high)),
                width = 0.4, size = 1,
                show.legend = FALSE)+
  ylab("Predicted Probability of Consent")+
  xlab("Recruitment Method")+
  scale_y_continuous(labels=scaleFUN)+
  scale_color_manual(name = "Recruitment path",
                     values = c("#000000", "#000000"))+
  ggh4x::facet_grid2(. ~ outcome,
                     scales = "free", independent = "y", switch = "y", remove_labels = "x")+
  ggh4x::scale_y_facet(COL == 1, limits = c(0.35,0.60), breaks = c(.4,.5,.6), labels = scales::percent) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none")

py <-
  ggplot(df_emm_installation_F4, aes(x = rev(as.factor(recruitment_arm)), y = rev(ame_predict), color = rev(as.factor(recruitment_arm)))) +
  geom_point(size = 5) +
  geom_errorbar(data=df_emm_installation_F4,aes(ymin=rev(conf.low),ymax=rev(conf.high)),
                width = 0.4, size = 1,
                show.legend = FALSE)+
  ylab("Predicted Probability of Installation")+
  xlab("Recruitment Method")+
  scale_y_continuous(labels=scaleFUN)+
  scale_color_manual(name = "Recruitment path",
                     values = c("#000000", "#000000"))+
  ggh4x::facet_grid2(. ~ outcome,
                     scales = "free", independent = "y", switch = "y", remove_labels = "x")+
  ggh4x::scale_y_facet(COL == 1, limits = c(0.35,0.60), breaks = c(.4,.5,.6), labels = scales::percent) +
  theme_bw(base_size = 24) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none")


pz <-
ggplot(df_emm_active_F4, aes(x = rev(as.factor(recruitment_arm)), y = rev(ame_predict), color = rev(as.factor(recruitment_arm)))) +
  geom_point(size = 5) +
  geom_errorbar(data=df_emm_active_F4,aes(ymin=rev(conf.low),ymax=rev(conf.high)),
                width = 0.4, size = 1,
                show.legend = FALSE)+
  ylab("Predicted Number of Tracked Days")+
  xlab("Recruitment Method")+
  scale_y_continuous(labels=scaleFUN)+
  scale_color_manual(name = "Recruitment path",
                     values = c("#000000", "#000000"))+
  ggh4x::facet_grid2(. ~ outcome,
                     scales = "free", independent = "y", switch = "y", remove_labels = "x")+
  ggh4x::scale_y_facet(COL == 1, limits = c(24,40)) +
  theme_bw(base_size = 10) +
  theme(legend.position = "none")


grid.arrange(px,py,pz, nrow = 1)

```







